#include <fstream>
#include "LoadTetFile.h"
#include <algorithm>
using namespace std;

bool fn_LoadGdmFile(std::string &sName, std::vector<Node> &vR, std::vector<Tria> &vElem)
{
    std::string sFile = sName + ".gdm";
    std::ifstream fp(sFile);
    if (!fp.is_open()) {
        return false;
    }
    int iNumNode, iNumWire, iNumTria, iNumQuad, iNumBoud;
    int m, n, k;
    fp >> iNumNode >> iNumWire >> iNumTria >> iNumQuad;
    fp >> iNumBoud;
    for (auto i = 0; i < iNumBoud; ++i) {
        fp >> m >> n >> k;
    }
    vR.resize(iNumNode);
    vElem.resize(iNumTria);
    for (auto &r : vR) {
        fp >> r.x >> r.y >> r.z;
    }
    for (auto &tri : vElem) {
        fp >> tri.node[0] >> tri.node[1] >> tri.node[2] >> tri.id;
        --tri.node[0];
        --tri.node[1];
        --tri.node[2];
    }
    fp.close();
    return true;
}

void fn_MarkGdmElem(std::string &sName, std::vector<Node> &vR, std::vector<Tria> &vElem, std::vector<bool> &vRefine)
{
    Node rc;
    Real dRadius, dist;

    vRefine.reserve(vElem.size());

    rc.x = 0.0;
    rc.y = 0.0;
    rc.z = 0.0;
    dRadius = 0.3;

    Node d;
    for (auto& tri : vElem) {
        d.x = (vR[tri.node[0]].x + vR[tri.node[1]].x + vR[tri.node[2]].x)/3.0;
        d.y = (vR[tri.node[0]].y + vR[tri.node[1]].y + vR[tri.node[2]].y)/3.0;
        d.z = (vR[tri.node[0]].z + vR[tri.node[1]].z + vR[tri.node[2]].z)/3.0;
        d.x -= rc.x;
        d.y -= rc.y;
        d.z -= rc.z;
        dist = sqrt(d.x * d.x + d.y * d.y + d.z * d.z);
        if (dist < dRadius && dist >0.2) {
            vRefine.push_back(true);
        } else {
            vRefine.push_back(false);
        }
    }
}

bool fn_OutputMesh(std::string &sName, std::vector<Node> &vR, std::vector<Tria> &vElem)
{
    std::string sFile = sName + ".msh";
    std::ofstream fp(sFile);
    if (!fp.is_open()) {
        return false;
    }
    fp << "MESH dimension 3 ElemType Triangle Nnode 3" << std::endl;
    fp << "Coordinates" << std::endl;
    for (auto i = 0; i < vR.size(); ++i) {
        fp << i + 1 << " " << vR[i].x << " " << vR[i].y << " " << vR[i].z << std::endl;
    }
    fp << "End Coordinates" << std::endl << std::endl;

    fp << "Elements" << std::endl;
    int v[4];
    for (auto i = 0; i < vElem.size(); ++i) {
        v[0] = i + 1;
        v[1] = vElem[i].node[0] + 1;
        v[2] = vElem[i].node[1] + 1;
        v[3] = vElem[i].node[2] + 1;
        fp << v[0] << " " << v[1] << " " << v[2] << " " << v[3] << std::endl;
    }
    fp << "End Elements" << std::endl;
    fp.close();
    return true;
}


bool fn_LoadTetFile(std::string &sName, std::vector<Node> &vR, std::vector<Tetr> &vElem)
{
    std::string sFile = sName + ".tet";
    std::ifstream fp(sFile);
    if (!fp.is_open()) {
        return false;
    }
    int iNumNode, iNumTetr, iNumBoud;
    fp >> iNumNode >> iNumTetr >> iNumBoud;
    vR.resize(iNumNode);
    vElem.resize(iNumTetr);
    for (auto &r : vR) {
        fp >> r.x >> r.y >> r.z;
    }
    for (auto &tet : vElem) {
        fp >> tet.node[0] >> tet.node[1] >> tet.node[2] >> tet.node[3] >> tet.id;
        --tet.node[0];
        --tet.node[1];
        --tet.node[2];
        --tet.node[3];
    }
    fp.close();
    return true;
}

void fn_MarkTetElem(std::string &sName, std::vector<Node> &vR, std::vector<Tetr> &vElem, std::vector<bool> &vRefine)
{
    Node e1, e2, norm;

    auto vArea = vector<double>(2, 0);
    vector<vector<double>> vvList;
    vvList.reserve(vElem.size());
    for (auto i = 0; i < vElem.size(); ++i) {
        e1.x = vR[vElem[i].node[1]].x - vR[vElem[i].node[0]].x;
        e1.y = vR[vElem[i].node[1]].y - vR[vElem[i].node[0]].y;
        e1.z = vR[vElem[i].node[1]].z - vR[vElem[i].node[0]].z;
        e2.x = vR[vElem[i].node[2]].x - vR[vElem[i].node[0]].x;
        e2.y = vR[vElem[i].node[2]].y - vR[vElem[i].node[0]].y;
        e2.z = vR[vElem[i].node[2]].z - vR[vElem[i].node[0]].z;
        norm.x = e1.y * e2.z - e1.z * e2.y;
        norm.y = e1.z * e2.x - e1.x * e2.z;
        norm.z = e1.x * e2.y - e1.y * e2.x;
        e1.x = vR[vElem[i].node[3]].x - vR[vElem[i].node[0]].x;
        e1.y = vR[vElem[i].node[3]].y - vR[vElem[i].node[0]].y;
        e1.z = vR[vElem[i].node[3]].z - vR[vElem[i].node[0]].z;
        vArea[0] = -fabs(norm.x * e1.x + norm.y * e1.y + norm.z * e1.z)/6.0;
        vArea[1] = i;
        vvList.push_back(vArea);
    }
    sort(vvList.begin(), vvList.end());

    vRefine.assign(vElem.size(), false);
    int iNum = vElem.size() / 10;
    int m;
    for (auto i = 0; i < iNum; ++i) {
        m = int(vvList[i][1]);
        vRefine[m] = true;
    }
}

bool fn_OutputMesh(std::string& sName, std::vector<Node>& vR, std::vector<Tetr>& vElem)
{
    std::string sFile = sName + ".msh";
    std::ofstream fp(sFile);
    if (!fp.is_open()) {
        return false;
    }
    fp << "MESH dimension 3 ElemType Tetrahedra Nnode 4" << std::endl;
    fp << "Coordinates" << std::endl;
    for (auto i = 0; i < vR.size(); ++i) {
        fp << i + 1 << " " << vR[i].x << " " << vR[i].y << " " << vR[i].z << std::endl;
    }
    fp << "End Coordinates" << std::endl << std::endl;

    fp << "Elements" << std::endl;
    int v[5];
    for (auto i = 0; i < vElem.size(); ++i) {
        v[0] = i + 1;
        v[1] = vElem[i].node[0] + 1;
        v[2] = vElem[i].node[1] + 1;
        v[3] = vElem[i].node[2] + 1;
        v[4] = vElem[i].node[3] + 1;
        fp << v[0] << " " << v[1] << " " << v[2] << " " << v[3] << " " << v[4] << std::endl;
    }
    fp << "End Elements" << std::endl;
    fp.close();
    return true;
}

bool fn_OutBoudMsh(std::string& sName, std::vector<Node>& vR, std::vector<Tetr>& vElem)
{
    constexpr int LocNode[4][3] = { {0,1,2}, {1,2,3}, {0,2,3}, {0,1,3} };

    if (vElem.empty()) return 0;

    int iNode, m, n, j, k;
    bool bFind;
    auto vvNode2Face = vector<vector<int>>(vR.size());
    // set first tetrahedron
    vector<vector<int>> vvLib;
    auto vTetFace = vector<vector<int>>(vElem.size(), vector<int>(4));
    auto vFace = vector<int>(3);
    for (int i = 0; i < 4; ++i) {
        vTetFace[0][i] = i;
        for (j = 0; j < 3; ++j) {
            iNode = vElem[0].node[LocNode[i][j]];
            vFace[j] = iNode;
            vvNode2Face[iNode].push_back(i);
        }
        sort(vFace.begin(), vFace.end());
        vvLib.push_back(vFace);
    }
    // check the other tetrahedron
    for (int i = 1; i < vElem.size(); ++i) {
        // check each face
        for (n = 0; n < 4; ++n) {
            vFace[0] = vElem[i].node[LocNode[n][0]];
            vFace[1] = vElem[i].node[LocNode[n][1]];
            vFace[2] = vElem[i].node[LocNode[n][2]];
            sort(vFace.begin(), vFace.end());
            bFind = false;
            for (j = 0; j < 3; ++j) {
                iNode = vFace[j];
                for (k = 0; k < vvNode2Face[iNode].size(); ++k) {
                    m = vvNode2Face[iNode][k];
                    if (vvLib[m][0] == vFace[0] && vvLib[m][1] == vFace[1] && vvLib[m][2] == vFace[2]) {
                        bFind = true;
                        vTetFace[i][n] = m;
                        break;
                    }
                }
                if (bFind) break;
            }
            if (!bFind) {
                m = int(vvLib.size());
                vTetFace[i][n] = m;
                vvNode2Face[vFace[0]].push_back(m);
                vvNode2Face[vFace[1]].push_back(m);
                vvNode2Face[vFace[2]].push_back(m);
                vvLib.push_back(vFace);
            }
        }
    }
    auto vFlag = vector<int>(vvLib.size(), 0);
    for (int i = 0; i < vTetFace.size(); ++i) {
        for (n = 0; n < 4; ++n) {
            ++vFlag[vTetFace[i][n]];
        }
    }

    std::string sFile = sName + "_boud.msh";
    std::ofstream fp(sFile);
    if (!fp.is_open()) {
        return false;
    }
    fp << "MESH dimension 3 ElemType Triangle Nnode 3" << std::endl;
    fp << "Coordinates" << std::endl;
    for (auto i = 0; i < vR.size(); ++i) {
        fp << i + 1 << " " << vR[i].x << " " << vR[i].y << " " << vR[i].z << std::endl;
    }
    fp << "End Coordinates" << std::endl << std::endl;

    fp << "Elements" << std::endl;
    int v[3];
    m = 0;
    for (auto i = 0; i < vvLib.size(); ++i) {
        if (vFlag[i] == 1) {
            ++m;
            v[0] = vvLib[i][0] + 1;
            v[1] = vvLib[i][1] + 1;
            v[2] = vvLib[i][2] + 1;
            fp << m << " " << v[0] << " " << v[1] << " " << v[2] << std::endl;
        }
    }
    fp << "End Elements" << std::endl;
    fp.close();
    return true;
}
